Mon. Not. R. Astron. Soc. 000, dlSl (2006) 



Printed 5 February 2008 



(MN WF^ style file v2.2) 



Turbulent diffusion in protoplanetary discs: 
The effect of an imposed magnetic field 



o 
o 

(N 



A. Johansen\ H. Klahr^ and A.J. Mee^ 



^ Max-Planck-Institut fiir Astronomie, Kdnigstuhl 17, 69117 Heidelberg, Germany 
^School of Mathematics and Statistics, University of Newcastle upon Tyne, NEl 7RU, UK 



Accepted 200- December -. Received 200- December — ; in original form 200- October ■ 



> 

in 

m 
o 

o 
o 



ABSTRACT 

We study the effect of an imposed vertical magnetic field on the turbulent mass dif- 
fusion properties of magnetorotational turbulence in protoplanetary discs. It is well- 
known that the effective viscosity generated by the turbulence depends strongly on 
the magnitude of such an external field. In this letter we show that the turbulent 
diffusion of the flow also grows, but that the diffusion coefficient does not rise with 
increasing vertical field as fast as the viscosity does. The vertical Schmidt number, i.e. 
the ratio between viscosity and vertical diffusion, can be close to 20 for high field mag- 
nitudes, whereas the radial Schmidt number is increased from below unity to around 
3.5. Our results may have consequences for the interpretation of observations of dust 
in protoplanetary discs and for chemical evolution modelling of these discs. 

Key v^rords: accretion, accretion discs — diffusion — MHD — turbulence — plane- 
tary systems: formation — planetary systems; protoplanetary discs 



1 INTRODUCTION 



Planets form out of micrometer-sized dust grains that 
are embedded in t he gas in protoplanetary discs (see 
for a recent review). The observed in- 
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frared radiation from protoplanetary discs comes primarily 
from micron-sized grains, although observations at longer 
wavelengths show that some discs have large populations of 
grain s with sizes up to mms and cms (e.g. Rodmanri et al. 
I2OO61) . Turbulent motions in the gas play a big role in the 
dynamics of chemical species and solids, at least as long as 
the solids are smaller than a few ten meters. Thus an under- 
standing of how dust grains and chemical species move un- 
der the influence of turbulence is vital for our understanding 
of the physical processes that take place in protoplanetary 
discs and the observational conseque nces ("llgn er et alJl2004l: 
l llgner fc NelsonlbOOd Iwillacv et alj |2006: Du Uemond et alJ 
boodlSemenov et aLlboO^ 

Turbulence has a number of effects on the embedded 
dust grains. Larger grains (rocks and boulders) can be 
trapped in the turbulent flow due to th eir marginal cou- 
pling to the gas (iBarge fc SommerialllQQSl) . whereas smaller 
grains feel the effect of the turbulence as a combination of 
diffusion and simple advection. Any bulk motion of the gas, 
e.g. turbulent motion with a turn-over time that is longer 
than the time-scale that is considered or even a radial accre- 
tion flow, leads to an advective transport of the dust grains 
rather than diffusion. The turbulent transport acts as diffu- 
sion only when the considered time-scale is longer than the 



eddy turn-over time. The turbulent diffusion coefficient of 
the grains, Dt = Sc^fi^^ , is often assumed to be equal to 
the turbulent viscosity of the gas ffow Vt = ac^fl^^ . Here 
a non-dimensionalisation with sound speed Cs and K eple- 
rian frequency Qq is used dShakura fc Sunvaevlll973) . The 
Schmidt number, a measure of the relative strength of turbu- 
lent viscosity and turbulent diffusion, is defined as the ratio 
Sc = ft/-Dt = a/ 5. Several recent works have measured the 
turbulent diffusion coefficient directly fro m numerical simu- 
lation s of magnetorotational turbulence iBalbus fc HawlevI 
I1991I) . The simulations by Johansen fc Klahr (2005, here- 
after JK05) yielded a Schmidt number that is around unity 
for radial diffusion, whereas Carballido, Stone, fc Pringle 
(2005, hereafter CSP05) found a value as high a s 10. The ver- 
tical S chmidt n umber, measured both by JK05,lTurner et al] 
J2006ri and by iFromang fc PapaloizoJ feOOd) . gives more 
consistently a number between 1 and 3. Here it is worthy 
of n ote that |Turiier_et_al . (2006) consider stratified discs, 
and iFromang fc Paoaloizou (.2006) even include the effect 
of a magne tically dead zone without turbulence a round the 
mid-plane jGammidllQiS iFleming fc Stonell2003h . 

This letter addresses the discrepancy between the dif- 
fusion properties of turbulence in protoplanetary discs re- 
ported in the literature. We show that a vertical imposed 
magnetic field affects the diffusion coefficient strongly. It 
is known that a net vertical field component leads to 
turbulence with a s tronger angular momentum transport 
iHawlev et al.|[^995^ . We perform computer simulations of 
magnetorotational turbulence for various values of the verti- 
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cal field and find that turbulent diffusion does not increase as 
much as the viscosity increases. Thus the ratio between vis- 
cous stress and diffusivity, i.e. the Schmidt number, also in- 
creases with the magnitude of the external field. As a result 
we are able to give a possible explanation for the discrepancy 
in the radial Schmidt numbers found in the literature. 



2 SOURCES OF AN EXTERNAL MAGNETIC 
FIELD 

The properties of any external magnetic field threading pro- 
toplanetary discs are not well-known. Close to the central 
object there is an interaction with the possibly dipolar or 
maybe quadrupolar magnetic field of the young stellar ob- 
ject. Also the occurrence of jet phenomena indicates that at 
least for the originating zone of the jet, e.g. a few protostar 
radii , there should be a large scale vertical magnet ic field 
re.g. iFendt Flstnetirmfli-lviemmings et al.ll2nnfj). How- 
ever, at larger orbital distances relevant for planet forma- 
tion, it is not obvious what the global field configuration 
should look like. 

To get some physical insight into the role of an external 
magnetic field in the dynamics of protoplanetary discs, we 
do here some rough estimations for two cases, either that 
the field originates in the central object, or that it comes 
from the molecular cloud core out of which the disc formed. 



2.1 Protostar 

The dipolar field of the central protostar dominates the 
gas pressure of the disc until a certain inner disc radius 
- Rin. This is typically a few times the protostel lar radius 
JCamenzindl Il990l : iKonigll Il99ll: IShu et alJ ll994D . Beyond 
J?in the interaction between the dipole field and the accre- 
tion disc is strongly unstable and le ads to an opening up 
of the pr otostellar dipole field lines (Mil ler fc Ston3 119971: 
iFendt fc Elstner 2000 : Kiikcr et al. 2003 ). Even if the pro- 
tostar could retain its dipolar field at larger orbital radii, 
the magnetic pressure exerted by the field lines would fall 
so quickly with orbital radius [Bl{r) oc r~^\ that it would 
be completely unimportant at several AU from the protostar 
where the gas planets are believed to form. 



2.2 Molecular cloud 

In molecular cloud co res the magnetic fi eld, Bcioud, can be 
as large as ~ 100 |J.G l* Bourke et al.ll200lh . The gas pressure 
in the disc can be written as P = c^p, where Cs is the sound 
speed and p is the gas density. The mid-plane density of an 
exponentially stratified disc with scale height H depends on 
the column density S as p = S /[\/2t^H). The scale-height 
to radius ratio H/r, which also corresponds to the ratio of 
local sound speed to Keplerian speed wk, can be used to 
rewrite the gas pressure at the mid-plane of the disc as. 
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The plasma beta of the external magnetic field is defined as 
the ratio between gas pressure and magnetic pressure 13 = 
f/-Pmag. One can write the following scaling for the plasma 



beta /3cioud due to the magnetic field from the molecular 
cloud, 



/3cioud = 5.9 • 10^ 



H/r\ ( mA ( Bciou, 



0.1 J \MqJ \ 



1 g cm" 



(lOOAu) • ''^^ 



Here /3cioud has a falling trend with r because the low gas 
density in the outer part of the disc makes the magnetic pres- 
sure more important there. For a sufficiently strong cloud 
field, the plasma beta could be relatively low at a disc ra- 
dius of several hundred astronomical units. 



3 SIMULATIONS 

We simulate a protoplan etary disc in the shearing 
sheet approximation ( e .g. ICo ldreich fc Trem aind Il978l : 
iBrandenburg et al.lll995l: iHawlev et alJll995^ . Here a local 
coordinate frame corotating with the disc with the Keple- 
rian rotation frequency Qq at a distance ro from the central 
source of gravity is considered. The coordinate system is 
oriented so that x points radially away from the central ob- 
ject, y points in the azimuthal direction parallel to the the 
Keplerian fiow, and z points normal to the disc along the Ke- 
plerian rotation vector fio- Numerical calculations are per- 
formed using the Pencil Code (a finite difference code that 
uses sixth order symmetric sp ace derivatives and a third or- 
der time-stepping scheme, see lBrandenbur3l2003D . 



3.1 Gas 

Considering the velocity field u relative to the Keplerian 
fiow u^'' — —{3/2)f2ox, the equation of motion of the gas is 



du 
'dt 



— + (u- V)-u - 



u:r' — = /(u) - V In p 



" dy 

+-J X {B + Boz) + f^ 
P 



(3) 



The left-hand-side of equation contains terms for both 
the advection by the velocity relative to the Keplerian fiow 
and for the advection by the Keplerian fiow itself. The terms 
on the right-hand side are the modified Coriolis force. 
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which takes into account that the Keplerian velocity pro- 
file is advected with any radial motion, the force due to the 
isothermal pressure gradient with a constant sound speed 
Cs, the Lorentz force (including the contribution from an im- 
posed vertical field of strength Bo) and the viscous force 
that is used to stabilise the numerical scheme. The viscosity 
term is a combination of sixth order hyperviscosity and a lo- 
calised shock capturing viscosity. The use of hyperviscosity, 
hyperdiffusion and hyperresistivity is explained in JK05. For 
the shock viscosity, where extra bulk viscosity is added in 
regions of flow convergence, we refer to lHaugen et alJ ll2004) 
for a detailed description. 

The evolution of the mass density is solved for in the 
continuity equation 
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Table 1. Measured turbulent viscosity and diffusion coefficients 



Run 


Lx Ly y. L z 


So 





a 




May 


Mas 


<5x 


Scx 






A 


1.32x1.32x1.32 


0.00 


oc 


0.0028 ± 0.0004 


0.053 


0.053 


0.041 


0.0031 


0.90 


0.0016 


1.75 


B 




0.01 


20000 


0.0078 ± 0.0015 


0.079 


0.092 


0.064 


0.0058 


1.34 


0.0031 


2.52 


C 




0.03 


2222 


0.0367 ± 0.0142 


0.197 


0.185 


0.140 


0.0225 


1.63 


0.0092 


3.99 


D 




0.05 


800 


0.1811 ± 0.0773 


0.416 


0.300 


0.181 


0.0574 


3.16 


0.0123 


14.72 


E 




0.07 


408 


0.5529 ± 0.0964 


0.761 


0.421 


0.330 


0.1984 


2.79 


0.0300 


18.43 


A4 


1.00x4.00x1.00 


0.00 


oo 


0.0015 ± 0.0002 


0.055 


0.036 


0.031 


0.0017 


0.88 


0.0009 


1.71 


B4 




0.01 


20000 


0.0038 ± 0.0009 


0.079 


0.057 


0.052 


0.0038 


1.00 


0.0024 


1.58 


C4 




0.03 


2222 


0.0414 ± 0.0176 


0.206 


0.182 


0.134 


0.0177 


2.34 


0.0078 


5.31 


D4 




0.05 


800 


0.0793 ± 0.0371 


0.279 


0.239 


0.179 


0.0268 


2.96 


0.0091 


8.71 


B4 




0.07 


408 


0.1242 ± 0.0694 


0.366 


0.291 


0.221 


0.0356 


3.49 


0.0121 


10.26 



|^ + u.Vp + 4°'|^ = -pV.K + /B, (5) 

where /d is a combination of sixth order hyperdiffusion and 
shock diffusion. The magnetic field evolves by the induction 
equation which we write in terms of the magnetic vector 
potential A, 

^+uf^^^^OoAyX + ux{B + Boz) + f^. (6) 

Again we use sixth order hyperresistivity and shock resistiv- 
ity, through the function in regions of strong flow con- 
vergence. The value of Bo sets the strength of an external 
vertical magnetic field. 



3.2 Dust particles 

The turbulent diffusion coefficient Dt of the flow is measured 
by letting dust grains settle to the mid-plane of the turbulent 
disc. The dust layer is represented as individual particles 
each with a position a;''' and velocity vector u^'' (measured 
relative to the Keplerian velocity u'y^y). The gas acts on 
a dust particle through a drag force that is proportional 
to but in the opposite direction of the difference between 
the velocity of the particle and the local gas velocity. The 
dust grains do not interact mutually and do not have any 
feedback on the gas. The equation of motion of the dust 
particles is 

= +5, (7) 

where the modified Coriolis force / is defined in equation 
|0J|, Tf is the friction time and g is an imposed gravitational 
field (see below). We assume in the following that ri is con- 
stant and thus independent of the relative velocity between 
the grain and the surrounding gas. In protoplanetary discs 
this is a valid assump tion for sufficiently small dust grains 
JWeidenschillinJl977H . We use a value of i7o7"f = 0.01 which 
is small enough that the diffusion coefficient should not dif- 
fer significantly from that of a passive scalar (which can 
be seen as a dust grain in the limit of a vanishingly small 
friction time) . This value is also large enough that the com- 
putational time-step is set by the Courant criterion for the 
gas and not by the friction force in the dust equations. 

The particles change positions according to the dynam- 
ical equation 



Under the effect of a special gravity field acting on the dust 
only, g in equation (|7|l, the particles fall either to the hori- 
zontal mid-plane of the disc, in the case of a vertical gravity 
field g — gz{z)z, or to a vertical "mid-plane" in the case of a 
radial gravity field g = gx{x)x. We use a sinusoidal expres- 
sion Qi = —go sin(fciii) with a wavelength that is equal to 
the size of the simulation box. In the equilibrium state, the 
sedimentation is balanced by the turbulent diffusion away 
from the mid-plane, and the dust number density n, for the 
case of a vertical gravity field, is given by (see JK05) 

lnn(2)=lnniH — cos(fczz) , (9) 

where ni is an integration constant. The equivalent expres- 
sion for the radial gravity case is found simply by replacing 
2; by a; in equation (|5J. 

We run simulations with different values of the external 
magnetic field strength Bo between and 0.07, correspond- 
ing to a /3 ranging from infinity down to approximately 400. 
Our computational unit of velocity is the constant sound 
speed Cs, length is in units of the disc scale-height H, and 
density is measured in units of mean gas density po . In these 
units the turbulent viscosity and the turbulent diffusion co- 
efficient, Ut and Dt, axe numerically equal to the dimension- 
less coefficients a and 5. The unit of the magnetic field is 
then [B] — Cs^/Jtopo and is chosen such that /xq = 1. For 
each value of Bo we run one simulation with a vertical and 
one simulation with a radial gravitational field on the dust 
particles. The diffusion coefficients Sx and 5z are found by fit- 
ting a cosine function to the logarithmic dust density. From 
the amplitude we then determine the diffusion coefficient 
using equation The run parameters and the results are 
shown in Table Two simulation box sizes are considered, 
a square box with a side length of 1.32 and an elongated 
b ox with (Lx, Ly, L z) — (1.0,4.0, 1.0) (similar to the setup 
of ISano et alJlioMl . For the first case we use a resolution 
of 64^ grid points and 1,000,000 dust particles. Simulations 
with 128^ grid points were done by JK05 and showed only 
small differences from the 64^ simulations in the measured 
Schmidt numbers. Each model is run for twenty local orbits, 
i.e. 20 X 27rJ7(^^, of the disc. The runs with an elongated box 
are done with 64 x 256 x 64 grid points and 4,000,000 dust 
particles. 



4 RESULTS 



(8) 



For each value of the imposed magnetic field we have mea- 
sured the Q- value from the Reynolds and Maxwell stress ten- 
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Figure 1. The Schmidt number plotted against the a value and 
the best power-law fit (dotted lines). The best fit has Sca; = 
4.6a°-26 and Sc^ = 25.3aO -*^. 



Figure 2. The correlation time of the turbulent mixing coeffi- 
cients versus the a-value. The correlation times fall significantly 
with increasing a. 



sors (see Table 0. The a-value grows approximately expo- 
nentially with Bo- An a-value close to unity can be reached 
already for Bq = 0.07 (corresponding to /3 ~ 400). A sim- 
ilar investigation into the depe ndence of a on an i mposed 
vertical field was undertaken bv lHawlev et al.l lll995l) . Com- 
paring with Table 1 in that work, one sees that there is a 
relatively good agreement between those results and ours. 
Magnetorotational instability with an imposed vertical field 
develops into a "channel" solution jHawlev fc Balbujll992t 
iGoodman fc Xulll994l: ISteinacker fc H enning"200l'l. charac- 
terised by the transfer of the most unstable MRI mode to the 
the largest scale of the simula tion box and the subse quent 
decay of this large scale mode (ISano fc Inutsukal200j) . Suf- 
ficiently strong vertical fields can even cause stratified discs 
to break up altogether llMiUer fc Stone"2000\ The creation 
and destruction of the unstable channel solution gives signif- 
icant temporal fluctuations in the measured stresses, evident 
in the standard dev iation of the t urbulent viscosity in Table 
□ (see also Fig. 1 of ISano fc Inut suka 2001). 

For measuring the turbulent diffusion coefficient we con- 
sider the logarithmic number density of the dust particles 
averaged from 10 to 20 orbits. We have chosen to calcu- 
late the diffusion coefficient directly from this average state, 
rather than calculating it from the instantaneous dust den- 
sity at a given time t, because large-scale advection fiow 
only works as diffusion when averaged over sufficiently long 
times. The average dust density was found to be in excel- 
lent agreement with the cosine distribution of equation @ 
with a deviation from a perfect cosine of less than 5% for all 
simulations. Thus diffusion is a good description of the tur- 
bulent transport over long time-scales. This is partly due to 
the fact that we consider diffusion at the largest scale of the 
ffow, i.e. at a scale that is similar to or larger than the en- 
ergy injection scale of the MRL Diffusion over length scales 
that are smaller than the energy injection scale should be 
weaker, because dust density concentrations at small length 
scales are not stretched by the full velocity amplitude of the 
larger scales, but only by the velocity difference that the 



larger scales exert over the much narrower dust concentra- 
tion. The exp(cos) equilibrium state for the dust density, 
however, has almost all of its power at the largest scale of 
the simulation box, so any scale-dependency of the diffusion 
coefficient should not have any influence on the equilibrium 
state (the fact that the logarithmic dust density in the equi- 
librium state is a cosine supports this). 

The measured turbulent diffusion coefficients are writ- 
ten in Table Q It is evident that the turbulent diffusion 
coefficient does not increase as fast with increasing verti- 
cal field as the turbulent viscosity does. In Fig. we plot 
the vertical and radial Schmidt numbers as a function of a. 
Both Schmidt numbers approximately follow a power law 
with a. Making a best-fit power law, we find the empirical 
connections 

Sc, = 4.6a"'% (10) 
Sc, = 25.3a°•*^ (11) 

Considering the two box sizes individually (black and grey 
symbols in Fig. 0, the radial Schmidt number is seen to 
rise slightly faster with increasing a in the case of the elon- 
gated box with Ly = 4, whereas the vertical Schmidt num- 
ber follows a trend that is independent of the box size. In 
ideal MHD simulations with /? = 400, CSP05 find a radial 
Schmidt number of around 10. Using a similar value for /3, 
we find that the radial Schmidt number rises from unity in 
the case of no external field to ~ 3 — 4 when (3 ~ 400. This 
may explain at least part of the discrepancy between the 
results by CSP05 and JK05. The box size used in CSP05 is 
1.0 X 6.28 X 1.0, and is thus comparable to our elongated 
box. We have tried with Ly = 6.28 as well, but found no 
significant difference in the results. 

It is interesting to note that iFromang fc PapaloizoiJ 
l)200ffl have an a-value of 0.015 and a vertical Schmidt 
number of 2.8. That fi t s alm ost perfectly in Fig. Since 
iFromang fc Paoaloizoij l)200(ll do not have an imposed ver- 
tical field in their simulations, this may mean that the rise 
in Schmidt number with a is something fundamental and 
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not only an effect of the imposed magnetic field, although 
further investigations would have to be done to explore this 
connection in more detail. 

4.1 Correlation times 

One can express the diffusion coefficient caused by the scale 
fc of a turbulent flow as Dk ~ u^ik- Here is the veloc- 
ity amplitude of that scale and £k is the typical length-scale 
over which a turbulent feature transports before dissolving. 
The advection length £k can be approximated by £k ~ u^tk, 
where tk is the correlation time, or life time, of a turbulent 
structure. Taking now an average (and weighted) correla- 
tion time Tcor of all the scales, one gets the mixing length 
expression for the diffusion coefficient in direction i, 

DI''> = Tcor^? , (12) 

valid for Fickian diffusion (for th e validity of Fickian diffu- 
sion see iBrandenburg et alj|2004) . Here the Mach number, 
-^/uf/ca, is the root-mean-square velocity fluctuation in real 
space. The diffusion coefficient should thus scale roughly 
with Mach number squared. We plot the correlation times, 
calculated from equation H12|l . of Sx and 5z versus the a- 
value of the ffow in Fig. |5| The correlation time of the tur- 
bulent diffusion coefficients falls steeply with increasing a- 
value, so even though the Mach number of the flow increases, 
the time a given turbulent structure has for transporting 
the dust becomes shorter and shorter. Since the correlation 
times of radial and vertical diffusion have approximately the 
same dependence on a, the ratio of the diffusion coefficients 
can be expressed as 5x/Sz ~ (Ma^/Ma^)^. The anisotropy 
in the diffusion coefficient in favour of the radial direction is 
then mostly an effect of the anisotropy between the radial 
and vertical Mach numbers. 



5 SUMMARY 

In this letter we report that the Schmidt number of magne- 
torotational turbulence depends strongly on the value of an 
imposed vertical magnetic field. For large values of the verti- 
cal field, the relative strength of the turbulent diffusion falls 
with respect to the turbulent viscosity. This could explain 
part of the discrepancy between measurements of the radial 
turbulent diffusio n coefficient in magnetorotational without 
an imposed field iJohansen fc Kla hr 2 005.) and with an im- 
posed field JCarbaUido et alJl2005ll . 
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